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MODELING COLLISIONAL CASCADES IN DEBRIS DISKS: THE NUMERICAL METHOD 
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ABSTRACT 

We develop a new numerical algorithm to model collisional cascades in debris disks. Because of 
the large dynamical range in particle masses, we solve the integro-differential equations describing 
erosive and catastrophic collisions in a particle-in-a-box approach, while treating the orbital dynamics 
of the particles in an approximate fashion. We employ a new scheme for describing erosive (cratering) 
collisions that yields a continuous set of outcomes as a function of colliding masses. We demonstrate 
the stability and convergence characteristics of our algorithm and compare it with other treatments. 
We show that incorporating the effects of erosive collisions results in a decay of the particle distribution 
that is significantly faster than with purely catastrophic collisions. 

Subject headings: methods: numerical ~ circumstellar matter - planetary systems - infrared: stars 



L INTRODUCTION 

More than 700 extrasolar planets have been identified 
to date in over 500 planetary systems Most of these 
planets were discovered via radial velocity measurements. 
As a result, only a handful of them are less than 10 
Earth masses; the vast majority are gas giants resem- 
bling Jupiter. They are also in extremely close orbits 
to their host stars, making these systems dramatically 
different from ours. A large number of additional can- 
didate transiting sy stems have been fou nd recently with 
the Kepler mission (|Borucki et al.|[20lol ). 

In contrast to the great majority of known exoplanet 
systems, our own solar system has a complex configura- 
tion with gas giants at significant distances from their 
central star and rocky planets/asteroids within the gi- 
ant planet zone. The direct detection of rocky planets 
and planctesimals around other stars is only feasible un- 
der very rare circumstances. One of the most produc- 
tive approaches is indirectly via the thermal emission of 
their planetary debris dust belts. Ever since the discover- 
les with IRAS (jA umann et al.lll984lBackman fc Parescd 
Il993f ). we know that extrasolar systems may harbor disks 
of dust/debris that are generated by planetesimal colli- 
sions. The dust reprocesses the stellar light and emits it 
as thermal radiation in the 10-1000 /.tm wavelength range. 
A prototypical example of such a system is Fomalhaut, 
where a planet is shepherding the star's debris disk re- 
solved in both sca ttered light dK alas et ar'2008") and in 
infrared emission ( Holland et alJ ' 2003; Stapclfeldt et al.l 
12004 iMarsh et al.fe oOSl. Debris disks highlight the con- 
stituents of planetary systems that are many to hundreds 
of AU away from their stars. 

With the launch of the Spitzer Space Telescope, many 
observations have been obtained to detect and possibly to 

^ also at Steward Observatory, University of Arizona, Tucson, 
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resolve debris disks in the infrared regime. Debris disks 
have been probed around all types of stars, both in stel- 
lar clusters and in the field. These observations showed 
that even though debris disks are common around stars 
of all spectral types, they are more likely to be detecte d 
in the earlier stages of stellar evolution (jWvattl [2008D . 
We have also learned that debris disks may be located 
close to or far from their central stars, that there are 
systems with multiple debris rings (such as our solar 
system), and that there can be wide v arieties of miner- 
alogi cal compositions within the disks (|Carpenter et al.l 
120091 ). Debris disk studies are now a major compo- 
nent of the Hers chcl observing program ()Matthews et al.l 
l2010t lEiroa et~ al. 2010), which will provide substantial 
advances in our understanding of their outer zones. 

Interpreting these results demands theoretical insights 
in a variety of areas. For example, attempts have been 
made to understand the evolution of debris disks as a 
function of stellar type by studying them i n stellar clus- 
ters of different ages. As concluded in iGaspar et al 

f09l), so lar- type stars in the field (iBeichman et al 
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clusters ( Gorlova et arTT200& amli iSiegler etall [2007[ ) 



may show a faster decay tr end compared to that ob- 
served for earlier- type stars (jRieke et al.l 120051 : iSu et al.l 
i2006), although the difference is subtle and needs con- 
firmation. The decay trends of the fractional lumi- 
nosity (fd = Lexc/L ^:) show a large range in values. 
iSpangler et"all (|200lD find a decay oc t"^ '^^ when fit- 
ting ISO/IRAS data, while IGreaves fc WvattI (|200l get 
a much shallower decay oc f~ 0-5 The majority of surveys 
however find a decay (x t~^ ()Liu et al.l 120041 : iMoor et all 
120061 : iRieke et all 120051 ). A better theoretical under- 
standing is needed to sort out these results and to pro- 
vide testable hypotheses that can be compared with the 
observations. 

Only a handful of debris disks have been resolved; for 
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the majority, we only know the integrated infrared excess 
emission. Finding the underlying spatial distribution of 
the debris in these disks is not straightforward, as any 
spectral energy distribution (SED) can be modeled with 
a degenerate set of debris rings at different distances. Al- 
though much of the uncertainty is associated with the op- 
tical constants of the grains, another under-appreciated 
issue is the grain size distribution. Collisional models can 
reduce the number of free parameters in the SED models 
by determining the stable size distribution of particles in 
the disks. 

Observations of resolved debris disks also have raised 
questions that can best be addressed by theoretical mod- 
els. For example, Spitzer MIPS images have shown a 
significant extended halo of dust around Vega ()Su et al.l 
[2)05), both at 24 and 70 /im. Initial calculations hypoth- 
esized the halo around Vega to be a result of a high out- 
flow of dust due to ra diation pressure from a recent high- 
mass c ollisional event ()Su et al.ll2005[) , while lMiiller et al] 
(|20T0t) model it as a result of weakly bound particles 
on highly eccentric orbits. Further modeling and deep 
observations of additional systems will help distinguish 
these two possibilities. 

In this paper, we describe a new algorithm for modeling 
debris disks, in which we refine the physics and numeri- 
cal methods used in collisional cascade models. In §2, we 
briefiy outline previous models and introduce the basics 
of our algorithm. In §3 we detail our numerical methods, 
followed in §4 by our approach for including simplified 
dynamics. In the last section we compare our numeri- 
cal algorithm to previous ones and discuss in detail the 
differences between the codes and the effects those differ- 
ences have on the outcome of the collisional cascades. We 
also supplement our paper with an extended appendix 
that covers the numerical methods and the verification 
tests of our code. 

2. THE PHYSICAL AND NUMERICAL CHALLENGES OF 
MODELING DEBRIS DISKS 

Collisional cascades have been studied both analyt- 
ically and using collisional integro-differenti al numeri- 
cal ni o dels. The clas sic an alytic models o f iDohnanvil 
(IT969I) . IHellved (Il970l ). and IBandermannI (I1972D took 
into account both erosive and catastrophic collisional 
outcomes, assumed a material strength that was inde- 
pendent of the particle mass, a particle mass distribu- 
tion with no cut-offs, and a constant interaction veloc- 
ity. These models yielded steady state power-law mass 
distribution indices of -11/6. This result was in gen- 
eral agreement with the measured size distribution of 
asteroids in the solar system . More rece ntly, analytic 
models by |Dominik fc De cipl ([200^ and iWvatt et al.l 
()2007f ) showed that the fractional infrared luminosity in 
a collision-dominated steady-state system decays follow- 
ing at~^ power-law . As introduced later in this section, 
iLohne et al.l (I2008D find a diff erent value for this decay 
timescale. IWvatt et al.l ()2007D also derived a maximum 
mass and fractional luminosity as a function of age and 
distance from the central star, which they then used to 
classify systems with possible recent transient events. 

However, numerical models are needed to expand on 
these results. In the particular case of our Solar System, 
sophisticated numerical models were developed t o track 
the evolution of the largest asteroids (Creenberg et al.l 



119781) . They have been further improved to reproduce 
the observed wavy structure in the size distribution 
at the very highest mass es (e.g., lO'Brien fc Creenberg] 
120051 iBottke et ai1l2005t ). These models yield power- 
law distributions that deviate from the classic solution 
of Dohnanyi (1969), with certain regions steeper than it 
and others shallower. Using a steeper or shallower dis- 
tribution and extrapolating it to dust sizes can result in 
substantial offsets in the number of particles and thus in 
the infrared emission originating from them for a given 
planetesimal mass. Conversely, the particle size distribu- 
tion affects the underlying disk mass calculated from the 
observed infrared emission. 

A complete numerical model of collisional cascades 
would follow outcomes from all types of collisions, in- 
clude a kinematic description of the system, incorporate 
coagulation below certain thresholds, and do all this with 
high numerical fidelity. Although such a model has not 
yet been built because of its complexity, there are a num- 
ber of approaches in the literature that model collisional 
cascades down to particles of micron size, each with dis- 
tinctive strengths and weaknesses. 

The collis ional code ACE ha s been used in many stud- 
ies ( Krivo v eTall [2?M l2005l l2006l 120081: ILohne et al.l 
12008; .Mliller et al.ll201Cil) . It follows the evolution of the 
particle size distributions as well as the spatial distri- 
bution of the dust in debris disks. The code initially 
only accounted for collisions result ing in catastrophic 
outcomes, while the latest version (jKrivov 
iMiiller et al.ll201Q) includes erosive (cratering) events as 
wel l. The coll i sional outcome prescriptions are based on 
the IDohnanvil ()1969D particle-in-a-box model, but with a 
more elaborate description of material strengths in col- 
lision outcomes as well as the radiation force blowout. 
A strength of the code is that it calculates the dynami- 
cal evolution of the systems, as well. Since following the 
dynamical evolution of a system makes large demands 
on computer memory space and CPU speed, the code 
can only model the size distribution with a low num- 
ber of mass grid points; it originally used a first order 
Euler Ordinary Differential Equation (ODE) solving al- 
gorithm, but has been modified to include a more pr ecise 
one (A. Krivov, pr iv. comm.). iKrivov et al.l ()2006[) and 
IMiiller et al.l ()2010f ) applied this algorithm to debris disks 
in general and to the specific example of the Vega sys- 
tem. They followed the orbital paths of fragme nts and 
placed special emphasis on radiation effects. Lohn e et al.l 
(|2008) modeled debris disk evolution around solar-type 
stars and found, both with analytic and numerical anal- 
ysis, that the majority of physical quantities, such as 
the mass and the infrared luminosity, decrease with time 
as t~^'^ to t~^''^. This is in contrast wi th the observed 
decay found by some obervati ons ()Liu et al.l 120041 : 
iMoor et al.l l2006 HRieke et al.l |2005D. However, th e pop- 
ulation synthesis verification tests in ILohne et al.l ()2008f ) 
yield good agreement with the latest Spitzer observa- 
tions. 

Thebauh et al. 1* 20031 l2007f ) study the evolution of ex- 
tended debris disks with a particle-in-a-box algorithm. 
They include both catastrophic and erosive collisions and 
employ resolution and numerical methods similar to the 
ones implemented in the ACE code. They model the ex- 
tended disk structure by dividing the disk into separate, 
but interacting rings. Their model does not include dy- 
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n amics. 

iCampo Bagatin et alj ()1994[ ) showed that a series of 
wave patterns is produced in the mass distribution of 
particles when a low-mass cutoff is enforced, such as in 
the case of a radiation force blowout limit. This signature 
is produced by the other de bris disk nume rical models as 
well fThcbau h et al.|[2003t IT hcbauh & Augereaul 120071: 
iKrivov et al.,i.2006l:lLohne et al..,200&.Wvatt et al.lboilD . 
However, the conditions under which these waves are pro- 
duced have not been completely analyzed. ^Wva tt et al.l 
(j^Oll) do show that the amplitude and wavelength of the 
waves is collisional velocity dependent. Such strong fea- 
tures in the particle size distribution are not observed in 
the dust collected within the solar sys t em. T he interplan- 
etary dust flux model of iGrun et al.l ()1985D . which used 
in situ satellite measurements of the micro-meteroid flux 
in the solar system, and the terrestrial particle flux mea- 
surements of the LDEF satellite ( Love fc Brownlec 1993,) 
only show a single peak at ^ lOO/im in the dust distri- 
bution. However, these measurements detected particles 
that were brought inward from the outer parts of the 
solar system via Poynting-Robertson drag and particles 
removed from the inner parts of the solar system via radi- 
ation force blowout. Their results are reflections of more 
than a single parent distribution and of multiple physical 
effects. 

The dynamical code of iKuchner fc Star Id (|20100 mod- 
els the evolution and 3D structure of the Kuiper belt, 
with a Monte Carlo algorithm and a simple treatment of 
particle collisions. Their models predict that grain-grain 
collisions are important even in a low density debris ring 
s uch as our Kuiper b elt. 

IKuchner fc StarkI (|20100 and iDuUemond fc DominikI 
(j2005f ) both emphasize the strong effects of fragmenta- 
tion in their models. Miiller_et al. (2010) point out that 
including erosive (cratering) events is necessary for their 
models to repr oduce the observed surface brig htness pro- 
files of Vega. iThebault fc Augereaul (|2007D also show 
that a complete collisional treatment will result in sig- 
nificant deviations from the classic power-law solution. 

Our goal is to set up a numerical model that places 
special emphasis on investigating these issues. Our new 
empirical description of collisional outcomes avoids dis- 
continuities between erosive and catastrophic collisions 
and thus enables a more stable and accurate calculation. 
We also solve the full scattering integral, thus ensur- 
ing mass conservation and the propagation of the largest 
remnants of collision outcomes. Finally, we use second 
order integration and fourth order ODE solving methods 
to improve the numerical accuracy. Below, we outline 
the physical and numerical techniques we will employ. 
In the Appendix, we present verification tests of these 
treatments. 

2.1. Collisional outcomes 

In collision theory, two types of outcomes are generally 
distinguished: catastrophic and erosive (the latter also 
known as cratering). For an erosive collision (EC), the 
collisional energy is relatively small, resulting in one big 
fragment whose mass is close to the original target mass. 
In the case of a catastrophic collision (CC), both colliding 
bodies are completely destroyed. 

We illustrate these outcomes in Figure [T] In the first 
panel, we plot the outcome distribution of an erosive col- 



lision, where there is a single large X fragment and a dis- 
tribution of dust at much lower masses. The X fragment 
is over half the mass of the original target mass M. The 
redistributed mass is equal to the cratered mass plus the 
projectile mass. The largest fragment in the distribution, 
Y, is arbitrarily set to be 20% of the cratered mass. In 
t j3.3l we elaborate on the validity of this arbitrary value. 

In the second panel of Figure [T] we plot the outcome 
at the boundary case between catastrophic and erosive 
collisions, where the single largest fragment X is exactly 
half of the original target mass M. The redistributed 
mass is equal to the other half of the target mass M plus 
the projectile mass. The largest fragment in the distri- 
bution, Y, is arbitrarily set to be 20% of the cratered 
mass here as well (10% of the target mass). 

Finally, in the third panel of Figure [11 we plot the 
outcome of a super-catastrophic collision, where the tar- 
get and projectile masses are equal. The mass of the 
single largest f ragme nt X is given by the relation of 
iFuiiwara et al.l ()1977[ ). The redistributed mass is equal 
to M — X plus the projectile mass. The largest fragment 
in the distribution, Y, is arbitrarily set to be at 0.5X. 

In reality, there is no strict b oundary between catas- 
trophic and erosive collisions (|Holsapple et al.l [2002). 
The outcomes between these two extreme scenarios 
should be continuous. In laboratory experiments, how- 
ever, it is easier to test the extreme outcomes. In our 
model, we use the laboratory experiments to describe 
the extreme solutions and connect them with simple in- 
terpolations throughout the parameter space. We re- 
vise the currently used models to include an X fragment 
for both erosive and catastrophic collisions as a separate 
new gain term. In this treatment, the placement of the 
X fragments is grid size independent, further improv- 
ing precision and guaranteeing the accurate downward 
propagation of these fragments. We are able to express 
the loss term in a much simpler form, including colli- 
sions from both regimes. Previous models only included 
a full loss term for catastrophic collisions and removed 
fractions of particles for erosive collisions. 

The slope of the power-l aw particle redi stribution has 
been studied extensively. iDohnanvil ()1969.' ) used a sin- 
gle power-law value from the largest mass to the small- 
est. Later experiments have shown that a double (or 
even a triple) po wer-law distribution i s a more likely out- 
come (see, e.g., iDavis fc RvanI 119901 ). This has led to 
the widespread use of a double power-law for the redis- 
tribution in numerous collision models. We conducted 
numerical tests that demonstrated that there is negli- 
gible difference in using a wide range of slopes with a 
single power-law. The fact that the outcome does not 
depend on the b imodality of the redistribution has also 
been proven by IThebault et al.l ()2003l ). Therefore, we 
have used the simplest method of redistributing the frag- 
mented particles with a single power-law slope from the 
second largest Y fragment downwards, scaled to conserve 
mass. As a nominal value, we will use -11/6 as the redis- 
tribution slope. This is close to the initial value of -1.8 
used by Dohnanyi (1969). 

2.2. Incorporating the complete redistribution integral 

The classic solution to the collis ional evolution of an 
asteroid system involves solving the lSmoluchowskH (|1916D 
integro-differential equation. This was first done by 
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Figure 1. Illustration of the possible outcome scenarios of collisions. In all collisions, a single largest X fragment is created as well as a 
power-law distribution of fragments with a largest mass of Y . 



iDohnanvil ()1969D . Because erosive collisions re move only 
a small part of the target mass in a collision, IDohnanvil 
([1969 ) expressed the erosive removal term in a differential 
form. This is not appropriate for our case. Our system 
has well defined boundaries; thus a continuity equation 
cannot be used. The locality of the collisional outcomes 
in phase space is not certain either. 

Therefore, to solve the Smoluchowski equation for the 
problem at hand, we need to solve the full scattering 
integral. This is complicated numerically, as the inte- 
grations must extend over the entire dynamical range of 
^ 40 orders of magnitude in mass. To be able to perform 
accurate integrations over such a large interval, we need 
to use a large number of grid points and sophisticated nu- 
merical methods. To achieve this in a reasonable time, 
we chose to drop the radial dependence of the various 
quantities and to solve the equations under a "particle- 
in-a-box" approximation. With this approach, we lose 
radial and velocity information but gain accuracy. 

2.3. The effect of radiation forces 

Poynting-Robertson drag can be an effective form of 
removing particles from the disk, so we include it in our 
model. However, the strongest and most dominating ra- 
diation effect is the removal of particles via radial ra- 
diation forces. These act on orbital timescales and can 
remove or place particles on extremely eccentric orbits. 
This gives rise to the challenge of incorporating a radial 
dependent removal term into a particle-in-a-box model 
that does not carry radial information. In §3.1.1. and 
§3.1.2, we discuss our approach for incorporating these 
radiation effects. 

Stellar wind drag is an important dust removal ef- 
fect for late type stars s u ch as in the case of AU Mic 
(|Augereau fc Beustll200l iStrubbe fc Chiand I2006D . an 
Ml spectral-type star. We concentrate on modeling de- 
bris disks in early- and solar-type systems, so we chose to 
neglect the effects of stellar wind drag. We do not take 
into account the Yarkovsky effect either, as it is small in 
high-density debris disks compared to the other radiation 
effects. 

3. THE COLLISIONAL MODEL 

We now discuss our collisional code (CODE-M - 
collisional Disk Evolution Model), which solves the sys- 
tem of integro-differential equations that describe the 
evolution of the number densities of particles of different 
masses. The code includes outcomes from erosive (cra- 
tering) collisions and catastrophic collisions and qualita- 
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Figure 2. Calculated values for the radiation-force parameter /3 
around stars of spectral-type AO, FO, GO, G5 and KG. The thin 
double-dashed black line is at the critical value /3=0.5, above which 
radiation forces are able to remove particles from circular orbits. 



tively follows the effects originating from radiation forces 
and Poynting-Robertson drag. 

Our system-dependent parameters are: the spectral- 
type of the central star (which defines the stellar mass 
M<t and the magnitude of the radiation effects it will 
have on the particles) , the minimum and maximum par- 
ticle masses (mmin and mmax, respectively), the radius, 
width, and height of the debris ring (i?, Ai?, and ft-, re- 
spectively), the total mass within the ring (A^tot), and 
the slope of the initial size distribution of the particles 
(77). We estimate the total volume of the narrow ring, V, 
as 

V = 2TThRAR, (1) 
which together with A^tot defines a mass density. 

3.1. The evolution equation 

In general, the change in the differential number den- 
sity n{m, t) at any given time for a particle of mass m is 
given by (j S moluchowskii [19161 ) 



dt 



n{m,t) = TpRD + T'coii , 



(2) 



where Tprd is the Poynting-Robertson drag (PRD) term 
and Tcou is the sum of the collisional terms. We define 
the differential number density of particles such that 



N{t) = / n{m,t)dm 



(3) 
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is the time-dependent total number density of particles 
within the ring. 

Effects such as radiation force blowout and Poynting- 
Robertson drag are able to deplete the low-mass end of 
the distribution, which in turn alters the evolution of 
the disk and more importantly, its infrared signature. 
Because we do not follow the radial profile of the vari- 
ous debris disk quantities in our algorithm, we can only 
capture the effects of radiation forces in a simplified way. 

3.1.1. Poynting-Rohertson drag term 

A complete analysis of the effects of the Poynting- 
Robertson drag is given bv I Burns et al.l (|1979t ). who cor- 
rect many errors made in previous work. This effect 
causes the parti cles to slow in th eir orbit and follow an 
inward spiral. Bu rns et al.l (|19791 show that the change 
in the orbital distance can be written as 



di?(m) _ 2GM,/3(m) 



cR 



(4) 



where G is the gravitational constant, c is the speed of 
light, and /3(m) is the ratio of radiation to gravitational 
force experienced by a particle of mass m. 

We calculate the f3{m) values as a function of the par- 
ticle masses, optical cons tants, and the spectr al type of 
the central star following iGaspar et al.l (|2008[ ). For the 
calculations we assume a silicate composition for the par- 
ticles and a bulk density of 2.7 g cm"'^. We show the 
calculated /3{m) values for a few different spectral type 
stars in Figured] 

We use equation ([¥]) to derive an approximate term 
that captures the effect of Poynting- Robertson drag as 
(see eq. [2]) 

_ n{m,t) 
JpRD — j-rr ' {^) 



TPRD {m) 



where 



TpRD (m) = 



2GM*/3(to 



■RAR 



(6) 



The mass dependence of the timescale comes from the 
mass dependence of the parameter /3. In principle, once 
a particle is removed from our coUisional system it still 
radiates in the IR; it just does not take part in the coUi- 
sional cascade. We keep track of the removal rate of these 
particles, but do not follow the total amount removed or 
their infrared emission. 

3.1.2. Radiation force blowout 

The effects of the radiation force blowout are incorpo- 
rated in our code with the simplified dynamics treatment 
introduced in §4, and not by the inclusion of a sepa- 
rate term in the differential equation as are the effects 
of Poynting-Robertson drag. Removing a particle from 
the coUisional system via radiation force blowout requires 
roughly an orbital timescale 



trfb = 27ri 



'_R3_ 
GM, 



(7) 



As we will show in §4, under our assumptions a newly 
created particle of mass m gets removed via radiation 
force blowout if /3(m) > 0.5 and is unaffected by radia- 
tion forces when /3(m) < 0.5. 
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Figure 3. Comparison of tiie radiation force blowout (RFB) to 
tlic Poynting-Robertson drag (PRD) timescales for an AO spectral- 
type star, as a function of particle mass and distance from the star, 
with a disk width of dR/R = 0.1. The dashed lines give the or- 
bital distances as a function of particle size, where the Poynting- 
Robertson drag and blowout timescales are 0.1, 0.01, 0.001 and 
0.0001 years. The solid red line gives the distance where the 
timescale for Poynting-Robertson drag is equal to the radiation 
force blowout timescale. Above the solid red line, radiation force 
blowout dominates, while below it Poynting-Robertson drag does. 
The plot shows that within reasonable disk radii estimates, radia- 
tion force blowout will always dominate in the li(rn) > 0.5 domain, 
while outside of it Poynting-Robertson drag will be the stronger 
effect. 



Although the radiation force blowout timescale is 
not used in our code, in Figure [3] we compare it to 
the Poynting-Robertson drag timescale around an AO 
spectral-type star, assuming a disk width-to-radius ratio 
of 0.1. The plot shows that within reasonable disk radii 
estimates, radiation force blowout will always dominate 
in the j3{m) > 0.5 domain, while outside of it Poynting- 
Robertson drag will be the stronger effect. Whether 
Poynting-Robertson drag is an effective form of removal 
in the /3(m) < 0.5 domain depends on the number den- 
sity of particle s in the ring, i.e., the coUisional timescale 
of the system (jWvattI I2005D . The outcomes are similar 
for all spectral type stars and for realistic AR/R values. 

3.2. The coUisional term 

The probability of a collision between particles is a 
function of their number densities and their coUisional 
cross section. We express the coUisional cross section for 
particles of mass m and m' as 



a (m, m') — t: [r (m) + r {m')] 

is the radius of each particle, 
rate of collisions between the 

P{m, m! ]t) = n(m, t)n{m! , t)V(7 (m, m!) 



(8) 

where r(m) is the radius of each particle. We express the 
differential rate of collisions between the two masses as 



; n(m, t)n{m! , t)V'n [r (m) + r (m')] 
: Kn{m, t)n{m' , t)VTT [m^ -\- m' 



(9) 



where V is their characteristic coUisional velocity, n{m, t) 
and n{m' , t) are the differential number densities for par- 
ticles of mass TO and to', 



(-Y 



(10) 
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Figure 4. The outcome possibilities as a function of colliding masses plotted for coUisional velocities of 0.5, 1 and 1.5 km s~^. These 
collisional velocities roughly correspond to debris ring radii of 100, 25 and 10 AU around an A spectral type star, respectively (see §4). We 
note that higher collisional velocities can occur in some systems. 



and p is the bulk mass density of the particles. The 
number densities in the problem are naturally all time 
dependent. However, for brevity, hereafter we drop the 
time dependence from our notation. 

The decrease or increase in number density at a cer- 
tain mass will be determined by three separate events: 
the removal of particles caused by their interaction with 
all other particles, the addition of the X particles from 
the interactions of other particles (see §2.1 and Figure 
[T]), and the addition of particles from the redistribution 
of smaller fragments originating from collisions of other 
particles. 

We express the first event, which describes the removal 
of particles, as 



-n{m) 



dm' Plm, m!) 



(11) 



We completely remove all particles from all grid points if 
they take part in a collision, even if they are the target 
objects in erosive collisions. 

The second event to be described is the addition of the 
large X fragments. To this end, we need to calculate the 
mass M that will produce a particle of mass m = X when 
interacting with a particle of mass m' . We achieve this 
with a root finding algorithm from the collisional equa- 
tions presented in the following sections and calculate it 
only once in the beginning of each run. In equation form. 



At 



n{m) 



tlx (m) 



dm'P{M, m') 



(12) 



m^X(m',Af) mmin 



The lower limit of the integration is the minimum mass 
in the distribution. We denote the largest mass m' that 
can create a particle of mass m as the X fragment as 
A*x(™)- Its value can also be calculated via root finding 
algorithms and has to be calculated for each value of m 
once in the beginning of each run (see Figure \W\ in the 
Appendix) . 

These first two integrals may catastrophically cancel, 
meaning that the difference between the two terms may 



be significantly smaller than the absolute value of each, 
causing the former to be artificially set to zero when eval- 
uated numerically. It is therefore useful to combine these 
terms into a single integral in a way that will lessen the 
probability of catastrophic cancellation: 



fix (ni) 



Ti{m) = — yTTK-j J dm'n(m') 



nim) 



l\ 2 

3 -I- m'^ 



-I- j drn'n{m')n{m)i^mi + m'^^ S , (13) 

fix{m) 

Unfortunately Ti can still suffer from catastrophic can- 
cellation (when m' is much smaller than m, and by defi- 
nition only for the first integral) . We overcome this issue 
by employing a Taylor-series expanded form of Ti, as 
given in the Appendix. 

The third event in the collisional term is the addition 
of the power-law fragments back to the distribution. The 
description of this process is quite simple; however, its 
precise calculation is not. We write in general 



Jii(m) 



d/i 

H[Y{fi,M) 



dMP{p, M)A{fi, M) X 



(14) 



where /i is the projectile mass, M is the target mass, A 
is the scaling of the power-law distribution, and H is the 
Heaviside function. The total redistributed mass is 



Mrcdist. (m, M) = 



Yifi.M) 



k{p,M)m-'<+^dm , (15) 



where Y is the largest fragment within the redistribution 
(i.e., the second largest fragment in the collision, after X; 
see §2.1). This gives the scaling factor 



(2-7)M,cdist.(M.M) 
r2-7 (^, M) 



(16) 
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The precision of this integration depends strongly on the 
resolution of our grid points, due to the integration limits 
set by the Heaviside function. We discuss in detail the 
integration methods we used in the Appendix. 

3.3. Collision outcomes 

The coUisional equations can be integrated if the val- 
ues of X{fi, M), Y{fi,M), and Mredist(/^, -M) are known 
as a function of the colliding masses. Their values are 
strongly dependent on the outcome of the collision they 
originate from, which is determined by the energies of the 
colliding parent bodies. We show the domains of erosive, 
interpolated erosive (explained later in the section), and 
catastrophic collisions as a function of the colliding body 
masses in Figure |4] for coUisional velocities of V — 0.5, 
1.0 and 1.5 km s~^. We introduce the method for calcu- 
lating coUisional velocities from orbital velocities in §4, 
but it is a good general approximation that the coUi- 
sional velocity is roughly an order of magnitude smaller 
than the orbital velocity. These coUisional velocities will 
then correspond to debris ring radii of 100, 25 and 10 
AU around an A spectral type star, respectively. 

A collision is considered to be catastrophic if 



2M 



> Q*{M) , 



(17) 



where Q*{M) is the dispersion strength parameter of the 
target mass M, fi is the projectile mass, and V is the 
relative velocity of the projectile compared to the par- 
ent rin g (§3.1). We use the disp ersion strength descrip- 
tion of iBenz fc Asphau^ (|1999[ ) and discuss our choice 
in the Appendix. Note that, in a more accurate treat- 
ment, we would redistribute the relative kinetic energy 
to both masses and not just to the target mass (i.e., di- 
vide hy fj, + M instead of M). We are, however, using 
the original definition of Qimpact (as opposed to using the 
relative kinetic energy) because the Q*{M) values that 
we will be com paring it to were defined the same way 
(jBenz fc Aspha ug 1999) and this definition makes the 
problem more tractable numerically. We note that some 
work has indicated that the tensile str ength curve itself 
may be collision velocity dependent (iBenz fc Asphaugj 
[19991: iHoisapple et al.ll200llStewart fc Leinhardtll2009D . 
which we currently do not take into account. 

In catastrophic collisions both particles are com- 
pletely destroy ed. Based on experimental evidence 
jFuiiwara et a ll 119771: IMatsui et al.l 119841: iTakagi et ahl 
[l984t ,Holsapple et al.ll200^ we will assume that apart 
from the largest fragment X{iJ,,M), the total mass is 
redistributed as a power-law distribution of particles 
up to a mass that we denote as F(/i, Af). We calcu- 
late the largest singl e mass created using the relation 
(jFuiiwara et al.lll97';n ) 



X{n,M) = M- 



2MQ*{M) 



-<9x 



(18) 



At the separatrix between catastrophic and erosive col- 
lisions, Q(/i, M)i^pact = Qd(A^)' and XipL.M) = M/2, 
which is exactly w hat we expect. The /3x factor is mea- 
sured to be 1.24 bv lFuiiwara et al.l ()1977D and this is the 
fiducial value that we use. Some experiments have shown 
that the shape and material of the target have an effect 



on th e exact value of /3x (jMatsui et al.lll98l : lTakagi et al.l 
I1984D . We will elaborate on the effects of varying /3x in 
an upcoming paper. 

The second largest fragment, Y , is always a fraction 
(0 < /y < 1.0) of the cratered mass, Mcr, in the erosive 
collision domain up to the erosive/catastrophic collision 
boundary. In the catastrophic collision domain, F is a 
fraction {fx) of the X fragment. We interpolate fx from 
its value defined by /y at the separatrix (where /y = 
fx, as X = Mcr) to a specified value Z™^" at the super 
catastrophic collision case of /i = M as 



fx = exp I ln(/y) 



In 



'max 
X 



In 



In 



2Q'{M)MV- 



M 

2Q*(M)My-^ 



(19) 

Our fiducial values for these fractions are /y = 0.2 and 
ymax _ express the remaining mass in catas- 

trophic collisions as 



Mrcdist(/^, M) ^ii + M- X(^, M) 



(20) 



which is redistributed as a large number of smaller par- 
ticles. 

Erosive collisions are more complicated and less well 
understood. A collision will be erosive if 



Q(/i,M) 



impact 



2M 



< QUM) 



(21) 



As described in §3, an erosive collision will result in a 
single large fragment, which will be a remnant of the 
target body, and a distribution of sm aller particles. We 
use the formula of iKoschnv fc GriinI (2001a,b,), i.e., 



(22) 



to calculate the total mass cratered from the target, 
where a and b are constants, with fiducial values of 
a = 2.7 X 10"^ and b = 1.23. This formula is only 
valid for small cratered masses; it can lead to artificially 
high values for the cratered masses (much larger than the 
target mass) even in the erosive collision domain. When 
the cratered mass given by this formula is larger than an 
arbitrarily set fraction /m of the target mass, we use the 
following interpolation formula 



Mcr = M X exp <^ ln(/M) + In 



0.5 \ 



fMj \n[QUM)IQi\ 



where 



a 



i/b 



(23) 
(24) 



We choose an arbitrary fiducial value for /m of 10"''. 

In erosive collisions, the single large fragment is ex- 
pressed as 

X = M - M„. . (25) 

As defined before, the largest fragment of the redis- 
tributed mass is a fraction /y of the cratered mass 



r(/X,Af) =/yA/er(M,M) , 



(26) 
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Table 1 

THE NUMERICAL, COLLISIONAL AND SYSTEM PARAMETERS USED IN OUR MODEL AND THEIR FIDUCIAL 

VALUES 



Variable Description 


Fiducial value Notes 


Numerical variable 


b Neighboring grid point mass ratio 


1.1 §C.2 


System variables 



p 

Mtot 

ri 

R 

AR 

h 

Sp 



Bulk density of particles 

Mass of the smallest particles in the system 

Mass of the largest particles in the system 

The total mass v^ithin the debris ring 

Initial power-law distribution of particle masses 

The distance of the debris ring from the star 

The width of the debris ring 

The height of the debris ring 

The spectral-typo of the star 

CoUisional variables 



§3, Eq. {TSll 
§3, Eq. lO 

§3 
§3 

Eqs. |[ni6l[7l [36ll40t 
Eqs. m l6l [361 l40t 
Eqs. lT1l37t 

S4 



7 Redistribution power-law 

$X Power exponent in X particle equation 

a Scaling constant in Mcr 

b Power-law exponent in Mcr equation 

/jvf Interpolation boundary for erosive collisions 

Jy Fraction of Y/Mcr 

ymax Largest fraction of Y/X at super catastrophic collision boundary . . . 

Constant in smoothing weight for large-mass coUisional probability . . 

P Exponent in smoothing weight for large-mass coUisional probability . 

Qsc The total scaling of the Q* strength curve 

S The scaling of the strength regime of the Q* strength curve 

G The scaling of the gravity regime of the Q* strength curve 

s The power exponent of the strength regime of the Q* strength curve 

g The power exponent of the gravity regime of the Q* strength curve . 



TTJW 

1.24 
2.7 X 10"'' 
1.23 
10-4 
0.2 
0.5 

lO^mmax 

16 

1 

3.5 X 10^ erg/g 
0.3 erg cm'^/g^ 
-0.38 
1.36 



Eqs 
Eqs 
Eqs 
Eqs 

Eqs. OS 
Eq. 
Eq. 1^ 



Eq. 
Eq. 
Eq. 
Eq. 
Eq. 
Eq. 
Eq. 



24 



Al 



while the redistributed mass is 



Mr 



(27) 



MrcdistifJ; M) ^ fl- 

Thus, the X{pL,M), Y{^,AI), and Mj-cdisti^J-, M) param- 
eters can be summarized as 

-/3x 



Mh 



2MQ' (M) 



M, 



rcdist 



[M ■ 

7x(M,M)X(Ai,M) 

JyAfer(/i,M) 

' fj. + M - X{p,M) 

^lJL + M„{lJL,M) 



in CC 
in EC 

in CC 
in EC 

in CC 
in EC 



(28) 



(29) 



(30) 



For the smallest particles, which we are particularly 
interested in modeling, radiation forces lead to effects 
such as reduced coUisional probabilities in thin ring disks 
and increased coUisional velocities in extended disks. In 
this section, we describe our approximate treatment of 
these effects. 

The radiation originating from the central star effec- 
tively modifies the mass of the star seen by the particles; 
the orbits themselves remain conic sections. The effec- 
tive mass of the star is decreased by a factor of 1 — /3(m), 
where /3(m) is defined in ij3.1.1l If 



(31) 



The Mrodist is redistributed as a large number of smaller 
particles in both collision types, with a slope of 7 = 
11/6 and a scaling given by equation ([T6|) . We choose a 
redistribution slope of 7 = 11/6, w hich is a value close to 
that given by experimental results (iDayis fc Rvanlll990D 
and is the same as used bv iDohnanvil ( 196^ 

We give a list of the variable coUisional parameters of 
our model and their fiducial values in Table [TJ 

3.4. The initial distribution and fiducial parameters 

We use the IDohnanvil ()1969[ ) steady-state solution of 
77 = 11/6 as our initial distribution, where 77 is the slope 
of the initial distribution and yields an initial number 
density of n{m) = Cm~^, where C is an appropriate 
scaling constant for the distribution. The exact value of 
this slope is unknown for all real systems. Fortunately, 
the convergent solutions and the timescales of reaching 
a convergent solution are fairly insensitive to this value. 

4. SIMPLIFIED DYNAMICS 



a newly created particle is generally put on a hyperbolic 
orbit, which we take a s the requirement for radiation 
force blowout to occur (|Kresaklll976l: IBurns et all 119791 
and references therein). 

The effects of dynamical evolution on the coUisional 
cascade can be traced to the eccentricity of the orbits. 
To follow the orbital path of a dust grain that has been 
created in a collision, we assume that the parent bodies 
were on circular orbits at radius R and had (3{m) ^ 0. 
We also assume that the produced grain will be created 
with very small relative velocity with respect to the par- 
ent bodies. Under these conditions, it can be shown that 
the semi-major axis of the acquired orbit will be 



, , l-/3(m) 
1 — 2p(m) 



(32) 



At /3(m) = 0.5 the semi-major axis becomes infinite, 
while at /3{m) = it is equal to the semi-major axis of 
the colliding particles' original orbit. The eccentricity 
(ep) of the orbit can be determined from the fact that 
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Figure 5. Left Panel: The weighted collisional probabilities of particles as a function of their mass m (and size). Right Panel: The weighted 
collisional probabilities of particles as a function of their eccentricity e. The particles are assumed to be within a narrow {dR/R = 0.1) ring 
at 25 AU from an AO spectral-type star. 



the periapsis will equal the original orbital distance 



yielding 



ep{m) 



a{m) [1 — efj{m)] — R 



P{m) / [I - P{m)] if;3<0.5 
> 1 if /3 > 0.5 



(33) 



(34) 



At /3(m) = 0.5, the eccentricity equals 1, and at /3(to) = 
0, it equals zero, consistent with our expectations. Simi- 
lar derivation s can be found elsewhere (e.g., Harwit 1963; 
lKresaUri976( ). A particle on an eccentric orbit will have 
a modified probability of interaction with other particles 
in the parent ring, which we address in §4.2. 

4.1. Collisional velocities 

iLissauer fc Steward (|1993[ ) give the velocity of a plan- 
etesimal relative to the other planetesimals in the swarm 
(i.e., the collisional velocity), averaged over an epicycle 
and over a vertical oscillation as 



V = Vorh 




(35) 



where e is the maximum eccentricity and i is the maxi- 
mum inclination in the system. This equation is valid for 
a swarm of particles in Rayleigh distributed equilibrium. 
This condition is true for a system in quasi-collisional 
equilibrium. We use this equation to estimate the colli- 
sional velocity of all particles, setting 



and 



'2R 
h_ 

2R ■ 



(36) 



(37) 



The smallest particles that are in highly eccentric or- 
bits will have varying velocities along their trajectories. 
However, when at their periapsis, they will have their 
original orbital velocities, as by definition they are on 
eccentric orbits due to their original periapsis velocity. 
Because of this, in our simplified dynamical treatment 
we only use a single collisional velocity for all particles, 
which is described by equation ([55]). 



4.2. 



Reduced collisional probabilities of /? critical 
particles 



Particles with /3(m) less than 0.5, but which are still 
non-zero, called /3 critical particles, are thought to pro- 
duce halos around debris disks via the highly eccentric 
orbits radiation forces p lace them on (|Thebault fc Wul 
l2008HMimer et al.llMof) . 

For a particle to go into an eccentric orbit, it must ac- 
quire a radial velocity component that is different than 
zero. In collisions, fragments will be ejected in all di- 
rections with a certain velocity distribution. Since the 
smallest fragments will tend to escape with the highest 
velocities (e.g., Jutzi et al. 2010), it is a fair question to 
ask whether thermalization of velocity vectors and their 
high values is a stronger effect in placing dust particles 
on higher eccentricity orbits compared to radiation ef- 
fects that reduce the effective stellar mass. 

To answer this question, we need to examine the ori- 
gin of the particles that contribute to the increase of the 
differential density in each mass grid point. We calculate 
T// only integrating in /i space, thus calculating the rate 
of increase of the differential number density of particles 
with mass m that originate from collisions with targets 
of mass M. Our calculations show that there is a pro- 
nounced peak in M roughly on the same scale (at most 
one order of magnitude higher) as m itself. That is, most 
particles originate from targets ~ 3 — 5x larger i n size 
than the particle itself. The results of lJutzi et al.l ()2010f ) 
clearly show that the velocities acquired by collision frag- 
ments at 1 /3 sizes are more than an order of magnitude 
lower than the collisional velocities, meaning that, in the 
most extreme case, a fragment will receive up to a few 
tenths of a km s^^ radial velocity compared to its 10- 
30 km s~^ orbital velocity. We can thus safely say that 
particles that are created with /3(m) values similar to 0.5 
tend to be placed on eccentric orbits by the radiation 
forces rather than being dispersed. These orbits will ex- 
tend out from the initial debris disk ring, preventing the 
particles from being destroyed and from them creating 
other particles. 

Our approach to calculating a weighting factor for each 
particle mass, determined by the fraction of its orbital 
period it spends in th e parent ring is similar to that of 
Thebault fc Wul ()2008D . The orbital time of a particle in 
an elliptical orbit as a funct ion of its distance from the 
center of mass is (|Tafilll985D 



cos 



aep 



{t - to) 



GM 



(38) 
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where to = is the initial time at periapsis, I is its dis- 
tance at time t from the center of mass, and we omit the 
m dependences of and a for clarity. We estimate the 
semi-major axis as 



R-AR/2 
1-ep 



(39) 



and we calculate the time At needed for a particle to 
reach the outer edge of the disk at I — R+dR/2. Dividing 
At by the half of the orbital period gives the weighting 
factor for each mass m as 



AR{2 - Cfi) ~ 2Rep 




ep [AR - 2R) 



AR [ep ~ 1) {AR-2Rep 



{AR - 2Ry 



(40) 



We plot these weighting factors as a function of the parti- 
cle mass and orbital eccentricity in Figure [5j When ana- 
lyzing the particle distributions, we only plot the number 
of particles within the parent ring, which we calculate as 



n^ingijn) — n{m)w{rn) 



(41) 



4.3. 



Reduced collisional probabilities of the largest 
particles 

The very last grid point in the domain of solution will 
only reduce its number density, as it cannot gain from 
larger masses either as an X fragment or from being part 
of a redistribution. The full phase space removal, intro- 
duced in i )3.2i causes its evolution time to become very 
small compared to all others and leads to a numerical 
instability. In order to avoid this, we multiply the colli- 
sional rates with a weight that smooths to zero for the 
largest particles 



cr^(m) 



1 — exp ( 



l-exp(-^) 



(42) 



for both the projectile and target particle. We chose Q to 
be a number a few orders of magnitude larger than mmax 
and use an arbitrary p = 16. The modified collisional 
rates, therefore, read 



P{m,m') — VTTKn{m)n{m')w{m)w{m') x 



13 + m 



(43) 



We discuss the implications of our choice of the weight 
function and of its parameters in §5.2. 

5. RESULTS 

As we discussed in §2, collisional cascades in debris 
disks have been studied extensively in the past decades, 
with many different analytic and numerical solutions to 
the problem. To demonstrate the similarities and differ- 
ences between our model and some earlier ones, we show 
in the following subsection the results of a few compar- 
ison tests. The system variables used by our code for 
these runs are summarized in Table [5J 

We compare our numerical model to three previ- 
ous well known algorithms, the particle-in-a-box code 
of PThcbault ct al.. (20031) . the dynamical code ACE 



Table 2 

PARAMETERS USED FOR COMPARISON 
MODELS 



Variable 


Comparison to 
Thebault (2003) 


Comparison to 
Lohne (2008) 


p (kg m-3) 


2700 


2500 


"Imin (kg) 


1.42x10^^1 


1.42x10^21 


"Imax (kg) 


1.78x101* 


4.20x101* 


Mtot (M®) 


0.0030221 


1.0 


V 


11/6 


1.87 


R (AU) 


5 


11.25 


AR (AU) 


1 


7.5 


h (AU) 


0.5 


3.4 


Sp 


AO 


G5 


PRD 


off 


off 



CKriv ov et ZI[200(l[2005tlLohne et al.][2(MlMiiller et al.l 

2010 ), and the ID steady-state solver code of 
IWvatt et all (|20ll . Although we do make an effort 
to model their systems as accurately as possible, a true 
benchmark between the codes is impossible. This is due 
in part to the fact that all models have somewhat differ- 
ent collisional and dynamical prescriptions. 

5.1. Comparison to \Thebault et al\ 1(20 OA ) 

A relatively straightfor ward comparison can b e made 
between CDDE-M and the iThebault et al.' ("20031) model. 
Although the initial IThebault et al. | (2003) model has 
been subsequently improved in IThebault fc Augereaul 
(|20f)7l ). we chose to compare our results with the former, 
as they are both particle-in-a-box approaches to the col- 
lisional cascade, with some dynamical effects included in 
a simplified mann e r. 

iThibault et al.l ()2003( ) aimed to model the inner 10 
AU region of the /3 Pictoris disk, with their reference 
model being a dense debris ring at 5 AU, with a width 
of 1 AU and height of 0.5 AU. We adopted their largest 
particle size of 54 km for the comparison run. How- 
ever, we adopted a smaller minimum particle size than 
they did (in our case well below the blowout size), to 
be able to follow the removed particles more completely. 
This choice does not affect the actual distribution within 
the ring. We also had Poynting-Robertson drag turned 
off. Although both of our models in clude erosive (cra- 
tering) collisions, the IThebault et all (2003) prescription 
uses hardness constants (a) of much softer material than 
that of our nominal case and a linear relationship be- 
tween cratered mass and impact energy (the prescrip- 
tion for erosive collisions has been changed in their later 
paper Thebault & Augereau l2007f ). For a better agree- 
ment, we also model a modified cratered prescription 
case, where we set b = 1, a = 10~* and fm = 0.01. 
With these adjustments our cratering prescriptions agree 
better; h owever our interpolati on formula is offset com- 
pared to [ThibauiFiFin dlOOl) . While ours has a con- 
tinuous prescription at the CC/EC boundar y (i.e., the 
cratered mass is 0.5M), the IThebault et al.l (|2003l) in- 
terpolatio n does not ( i .e., th e cratered mass is O.IM). 
Thebault fc Augereaul (|2007l ) improve on this, by em- 
ploying an interpolation formula very similar to our equa- 
tion ([jg]) . 

Figure [5] compares the evolution of th e distr ibution 
of particles between the IThebault et al.l (|2003[) nomi- 
nal case and our runs. In the vertical axes we plot 
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Figure 6. Comparison of tlie evolution of tlie dust distribution around the fi Pictoris disk modeled bylThebault et all (2003) and the model 
presented in this paper. The thin solid line is the initial distributi on. The modified cratered mass (Mcr) model uses an erosive collision 
prescription tha t is a closer analog to t he original IThebault et al.l l(2)03) soft material one. We also plot the curve from the innermost 
annulus of Thebault fc Au gereaul 1120071 ) as reference, scaled to the same density at the largest sizes as our model. We emphasize though 
that the Thebault fc Augereaur il2007') model is of an extended system, which is significantly different than the IThebault et al.l (12003) model 
or ours and is only plotted as a reference. It is noteworthy, however, that the dynamics and extendedness of the disk structure seemingly 
have less effect on the results than the material properties or the coUisional prescription itself. See text for more details. 
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Figure 7. Comparison of the evolution of the total disk mass and the dust-to-planetesimal mass ratio around the fi Pictoris disk modeled 

by Thebault et al. (2003) and the models presented in t his pap er. The modified cratered mass (Afcr) model uses an erosive collision 

prescription that is a closer analog to the original iThebault et al.l (120031 ) one. See text for more details. 



n(m) X m^, which is s i milar to the "mass/bin" value 
used bv IThebault et al.1 (I2003D. T o make them exact, we 
divide the IThebault et al.l (|200"3l) values by (5-1), which 
places them on the same scale. A few similarities and 
a few significant differences can be noted. Generally 



both models show wavy str ucture - which is a well stud- 
ied phenomenon (see e.g., ICampo Bagatin et al.l 1199^ 
IWvatt et al.ll201ll ) - but the exact structure of the waves 
differs. 

Our modified erosive (cratering) prescription model 
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Figure 8. Compar i son of ttie evolution of the particle distribution within a debris disk around a solar- type star modeled bv lLohne et al.l 
1)20081 ) . FWvatt et al.l 1)20111 ) and the code presented in this paper. The "only CC" CODE-M model uses only catastrophic collisions. The thin 
solid line is the initial distribution. See text for more details. 
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Figure 9. Comparison of the evolution of the dust mass within a debris disk around a solar-type star modeled bv lLohne et al.l 1)20081 ) and 
the model presented in this paper. See text for details. 



gives a much better agreement with the iThebault et alJ 
(|2003f ) results than our fiducial prescription, in the sense 
that it yields a deeper first wave in the distribution with 
a larger wavelength. The offset between the locations of 
the first dip and the subsequent peak in the two mod- 
els coul d likely result from the higher collisional veloci- 
ties that IThebault et al.l (|2003| ) calculate for the smallest 
particles. Our mo dified erosion prescription gives good 
agreement with the ,Thebault et a l. (2003) results for par- 
ticl es larger than a km in size, which is a surprise as 
the IThebault et al.l (|200l erosive constants are for much 
softer materials than our nominal values. Just above the 
blow-out regime our model becomes abundant in dust 
particles, as more and more dust is placed on highly ec- 
centric orbits. Although some smoothing is expected in 
reality, we do expect the number of dust particles near 
the blowout limit to increase. 

While both our distributions show the typical dou- 
ble power- law fea ture of quasi stead y-state collisional 
cascades (see e.g.. iWvatt et al.ll201lD above and below 
the change in t he str ength curve, it is masked in the 
IThebault et al.l ()2003[ ) model, due to the high ampli- 
tude wavy structures. These structures are substantially 
reduced for the innermost ring of the disk modeled in 
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time (yr) 

Figure 10. Evolution of the total mass within consecutive mass 
regions from the smallest to the largest particles in the system 
for the full collisional system, using the ii-0 . 3 parameters of 
ILohne et al.. t200a). The plot can be compared to the top panel of 
Figure 4. of ILohne et al.l 1)2001 ). 

IThebault &: Augereatj (|2007[ ). but persist in the outer 
zones. 
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Figure 11. Comparison of the evolution of our model distribution using the lLohne et al] 1)20081 ) ii-0.3 parameters, when varying the 
parameters of the weights in the collision cross sections for large particles. 



In Figure [71 we show the differences in the evolu- 
tion of the di s k ma ss between the nominal case of 
iThebault et"all (pOOl and our models. The left panel 
shows the evolution of the total disk mass within the de- 
bris ring, while the right panel shows the evolution of 
the dust-to-planetesim al mass ratio . These figures are 
equivalent to Thcbaul t et al.l (|2003D fi gures 2 and 3 (ex- 
cept that these are in plotted in logarithmic scales). Our 
nominal model predicts a faster decay of the total disk 
mass, reaching 25% mass loss, while the modified ero- 
sive prescription agrees with the Thebault et al. (2003) 
model and loses ~ 12% of its initial mass. The evolution 
of the dust-to-planetesimal disk mass differs while the 
quasi steady-state is being reached, after which all mod- 
els decay with the same slope. Our nominal case model 
has an order of magnitude larger dust-to-planetesimal 
mass ratio at all times compared to the pThebaul t et al.l 
(Hods') model, while our modified erosive collision pre- 
scription case is close to it. 

There are so me easily identi fied differences between 
our models. IThebault etahl (^20.03) use the same 
IBenz fc Asphaud (|1999D dispersive strength curve as we 
do, although they do average it to account for impact 
angle variations. This is an unnecessary step, as the 
IBenz fc Asphaud (|1999D strength curve is already impact 
angle averaged, and is corrected in lThebault fc Augereaul 
(|2007f) . However, we find that this scaling offset does 



not have a significa nt effect on the outcome of the dis- 
tribution evolution. IThebault et al.l ([2003) use a double 
power-law for fragment redistribution, while we use only 
a single power-law. We find that varying the slope of the 
single power-law does not have a significant effect on the 
evolution of the distribution either, so it appears that 
this difference is also likely not a significant contribut- 
ing factor to the discrepanci es. A noteworthy diffe rence 
between our models is that IThebault et al.l (|2003[ ) cal- 
culate fragment re-accumulation, while we do not. This 
is a possible explanation for our discrepancies at high 
masses, and the offsets we have in the total mass decay. 

The most significant difference between the models 
is that ours u s es a single interaction velocity, while 
IThebault et al.l (|2003D model the interaction velocity be- 
tween the /3 critical elliptical orbit smallest particles and 
the parent ring. This is likely to account for some 
of the additional offsets for the smallest particles, as 
higher interaction velocities have been shown to initi- 
ate higher amplitu d e waves (iCampo Bagatin et al.lll994l : 
IWvatt et al.ll20Tll) . IThebault et al.l (I2003D also take into 
account the constant presence of particles smaller than 
the blowout limit within the parent ring, which we do 
not. Within denser disks this step may smooth out any 
features near the blowout limit. 
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5.2. Comvaris on to \ Ldhne et all K200^ ) and 
VWvatt et al\ K2011 ) 

As introduced in §2^ th e numerical code ACE 

(IKrivov et al.l[2000l [200l [20061 [20081 ILohne et al.l[2008l: 
[Miiller et al.ll20ld[ ) solves the dynamical evolution of the 
coUisional system as well as the collisional fragmenta- 
tion, thus a straight comparison to CODE-M can not be 
performed. We use their ii-0.3 model ( Lohne et al.l 
[2008.) for comparison, which is for a relatively wide (7.5- 
15 AU), extremely dense (1 total mass with a largest 
planetesimal size of 74 km) debris ring. We turned off 
the effects of Poynting-Robertson drag for this compari- 
son model. The initial parameters we assumed are sum- 
marized in Table [2] This same system was modeled by 
IWvatt et'all (|2011[ ). whose results we also use for com- 
parison. 

In Figure [51 we show the evolution of the dust dis- 
tribution of the s ystem given by CODE-M, ACE, and 
Wvatt et all j2m^. As the version of ACE used in 



Lohne et al. (2008) o nly m odeled catastrophic collisions 



and the Wyatt et al] (|2011[ ) model uses catastrophic col- 
lision rates, we also include a CODE-M run in the plot 
that only models the outcome s of catastrophic collisions. 
Since the ILohne et al.l (|2OO80 values are already down- 
scaled by {6 — 1 ), no additional sca ling was required of 
their data. The IWvatt et al.l ()2011| ) data points are di- 
vided by (5 — 1 = 0.0626 to convert them to differential 
number densities. Qual itatively, all dis t ributi ons agree 
much better than in the iThebault et al.l (|2OO30 compar- 
ison (Figure [6]); however, there is some scaling offset be- 
tween the full CGDE-M and the other models, especially 
at large masses. 

The wavelengths of the waves roughly agree between 
CODE-M and ACE, with the single difference being the 
absence of the strong offset of the first crest in our 
model; the agreem ent is also good between CODE-M and 
IWvatt et al.] poTTI ). The double power-law distribution 
due to the change from strengt h to gravity dominated 
thresholds in the strength curve (jBenz fc Asphau3ll999| ) 
can be distinguished in all three models, with roughly 
the same slopes. The ACE and the IWvatt et all (|2011[ ) 
models maintain their initial —1.87 number density dis- 
tribution slope, while CODE-M becomes somewhat steeper 
for the smallest particles0 Our catastrophic-collision- 
only model has a smaller amplitude and wavelength wave 
structure than our full model or the other models. The 
most significant difference between the models is the scal- 
ing offset of the full CODE-M model, which we analyze 
below. 

In Figure [9l we sho w a comparison to figures 1 and 2 
of lLohne et al.l (|2008| ) . In the left panel we show the evo- 
lution of Cr, (CT,nVir i fi)i wh ich is introduced in Equation 
(11) of ILohne et all (|2008l) as 



Cr 



M. 



disk 



ML,. 



(44) 



This quantity is inversely proportional to the character- 
istic timescale of the system. As expected, since our sys- 
tem evolves faster, its characteristic timescale is shorter, 
so the Cl factor for our models is larger. This can be 

^ In Figure[8]we plot in the y-axis the product n(m)m^, so that a 
steeper number density slope will show up as a flatter distribution. 



seen in the right panel as well, where we plot the decay 
of t he total mass in ou r system and that given in figure 
2 of lLohiie et al.l (120081). We adopted the exact strength 
curve of ILohne et al.l (|2008[ ) in this run of our model, 
with the corrections given bv lWvatt et aD (poTTI ). 

In Figure [TOl we show a coni parison to the top panel 
of figure 4 in ILohne et al.l (|2008D , which shows the evolu- 
tion of the total mass within each of their mass bins. In 
this plot we show the evolution of the full collisional sys- 
tem, which includes erosive and catastrophic collisions. 
Since we do not use mass bins, but rather a differential 
number density griding, we integrate our distribution be- 
tween 14 grid points for each mass value, which roughly 
corresponds to a single mass bin of lLohne et all ([2008). 
Up to roughly a few hundred meters in size (where the 
strength curve has its minimum) all mass "bins" decay 
in close parallel slopes to each other after reaching their 
quasi steady-state around 10,000 yr. This is in contrast 
to the ILohne et al.l (I2008D res ults and agrees more with 
figure 2 of IWvatt et"ani20Tl[) . who model the same sys- 
tem. Our intermediate size planetesimals ( ~ km) show a 
steeper decay than that modeled by either ILohne et al.l 
(2008) or Wyatt ct al. (2011). 

The obvious difference between ACE and CODE-M, is that 
ACE also evolves the dynamics in the system and takes 
into account the varying collisional velocities in the sys- 
tem from particles that are in elliptical orbits within the 
parent ring. This could easily explain the offset of the 
first wave given by ACE in Figure [51 

The increasing offset between the full CODE-M run 
and the other two mod els is likely due to the omission 
of erosive collisions by ILohne et al.l (I2008D and to us - 
ing catastrophic colli s ion r ates in IWvatt eFaH ([20T1I) . 
iKobavashi fc Tanakal (|2010( ) have shown earlier erosive 
(cratering) collisions to be the dominant effect for mass 
loss in collisionally evolving systems. This effect is 
demonstrated by the CODE-M model we run with only 
catastrophic collisi ons included , which scales exactly 
with the ACE and IWvatt et all (|2011[ ) models. Since 
CODE-M does not include aggregation, the collisions of the 
smallest particles with the largest bodies is not modeled 
perfectly. We assume the realistic distribution decay to 
lie between the two models given by CODE-M. 

As introduced in ij4.31 we artificially reduce the colli- 
sion cross section of the largest particles in the system to 
zero in order to avoid numerical instabilities. However, 
as this is a completely arbitrarily defined numerical ne- 
cessity, we investigate its effects on the total mass decay, 
where we expect it to be the strongest. We reproduce 
Figures [51 & [TH but with varying the values of O and 
p in the smoothing function, in Figures [11] & [121 As 
can be seen in these plots, the variable O does not af- 
fect either the evolution of the distribution or the total 
mass decay, as long as it is larger than one. Varying the 
values of p does have an effect on both the evolution of 
the distribution and the total mass decay. The effect is 
only on the largest bodies in the system; below a size of 
one hundred meters, the shapes of the distributions re- 
main unchanged, with the only differences being scaling 
offsets. 

Visual examination of the distributions in Figure [HI hint 
at a slightly stee per distribut i on slo pe for the CODE-M 
models than the ILohne et al.l ([2008') and IWvatt et al.l 
(|2011i) ones. Since the distributions have wavy struc- 
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Figure 12. The change in the dust mass evolution when varying the parameters of the weights in the collision cross sections for large 
particles. 
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Figure 13. The evolution of the dust-to-pla netesimal ma s s ratio 
of the CODE-M models and the values for the IL5hne et al.l H2008I ') 
and IWvatt et al. (2011) models at 0.5 and 50 Myr. We calculate 
the dust mass from 0.1 mm to 10 cm and the planetesimal mass 
from 10 cm to 100 m. On the right vertical axis we give the value 
of the slope of the power-law distribution that would give the same 
dust-to-planetesimal mass ratio. 



tures in them, this is difficult to show with a slope fit. 
Therefore, we calculate the dust-to-planetesimal mass ra- 
tios of the distributions. We define the dust sizes to be 
from 0.1 mm to 10 cm and the planetesimal sizes to be 
from 10 cm to 100 m. In Figure [13] we plot the evolution 
of t he mass ratios of ou r mo dels and the niass r atios for 
the iLohne et all ()2008D and IWvatt et all (|2011[ ) models 
at 0.5 and 50 Myr. Our models have significantly higher 
mass ratios than theirs. This is likely a result of the 
differences between our coUisional equations. 



6. CONCLUSIONS 

In this paper wc present a numerical model of the evo- 
lution of the distribution of dust in dense debris disks. 
We calculate our model with a new numerical code, 
CODE-M, which we extensively verify and test in the Ap- 
pendix. 

Our coUisional model and numerical solution both 
present improvements to certain aspects of previ- 
ous numerical coUisional cascade models, but also 
have some limitations. Like most debris disk mod- 
els dXh cbauh & Augcrcau 2001 
Lohne et all 120081 iMiiUer et al] 

20111). we solve the Smoluchowski coUisional equation 



Krivov et al.l [20061: 
20101 IWvatt et al.l 



( Smoluchowskil [iQlGf) . which does not enable the study 
of stochastic impact events. Our model also does not 
evolve the dynamical sta te of the systems; the only code 
to currently do so is ACE (jKrivov et al.ll2005t ILohne et al.l 
[2008; Miillcr ct al. 2010). In its current state, our model 
is also limited to solving the coUisional cascade in de- 
bris rings, simi larly to the initial versions of most coUi- 
sional model s (iKrivov et al.|[2000l: iThebault et"all 120031: 
IWvatt et all 120111 ). some of which have since been ex- 
panded to model extended disk st ructures (jKrivov et al.l 
[2005t IThebault fc Au gerea ^ [20071 ). Ahhough our model 
only uses a single interaction velocity, as we show in H4.ll 
this is a valid assumption as long as the interactions oc- 
cur within the debris ring. 

We introduce a new prescription for describing erosive 
collisions, which always takes into account the reduction 
of the mass in the largest fragment during the cratering 
event. We introduce a number of approximate interpola- 
tions to ensure that our description of erosive and catas- 
trophic collisions yields a continuous set of outcomes as 
a function of the colliding masses, while being consistent 
with experimental results at various limits. Moreover, 
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we employ an efficient numerical algorithm that allows 
us to evaluate the scattering integrals with high preci- 
sion, considering the enormous dynamic range of masses 
involved in the calculation. 

We compare our code to the previously published nu- 
merical models of collisional cascades in debris disks, re- 
assuringly showing gene ral agreement , parti cularly with 
iLohne et all (|2008l ) and iWvatt et alj (|201lD . Nonethe- 
less, our mode l shows faster decays than previously 
published on es jThebault et alJ[200a ITohne et al.ll2008l : 
IWvatt et al.ll2011[ ) and also vields slightly higher dust-to- 
planetesimal mass ratios. We attribute these characteris- 



tics to be a result of our accurate treatment of collisional 
cascades. Our model will be used in a series of upcoming 
papers to study those aspects of debris disk behavior for 
which it is uniquely well suited. 

We thank Richard Greenberg and David O'Brien for 
helpful discussions on collisional outcomes and William 
Hartmann for advice on the smallest particles produced 
in collisions. We thank Philippe Thebault, Alexander 
Krivov and Mark Wyatt for their comments on an earlier 
version of the manuscript. Support for this work was 
provided by NASA through Contract Number 1255094 
issued by JPL/Caltech. 



APPENDIX 



A. STRENGTH CURVES 

The redistribution outcome of collisions depends almost solely on the energy of the impact and the colliding masses. 
In experiments it is common to specify the ratio of the k inetic energy of the projec tile to the mass of the target. This 
ratio is known as the specific energy Qimp of the impact. iGault fc Wedekiiidl ()1969t ) already noticed that the fragment 
distribution of particles depends on Qimp (which they called "rupture energy") when firing aluminum projectiles 
into glass targets. Their experiments showed that the fragments will have a power-law distribution, with the largest 
fragment being a fun ction of the specific energy of the impact. This r elationship was first given in equation format in 
iFuiiwara et al.l ()1977D for basalt targets. They note an offset from the IGault fc Wedekindl () 19690 results, likely due to 
material strength differences. 

Two specific values of Qimp are used: Q*g (the shattering specific energy) and a somewhat larger (the dispersion 
specific energy). The value of Q*g gives the energy required to shatter the target so that the mass of the largest 
fragment is no more than half of the original target mass. However, if the target is large enough, then self gravity 
pulls the fragments back together, leaving a remnant larger than half of the original. The larger gives the value 
of Qimp needed to disperse the fragments, so that the largest remaining piece is half of the original target mass. At 
lower target masses, where self-gravity can be neglected, QJi ~ Q*s- We use in our code and refer to it as Q* . 

Determining the value of Q* is difficult, especially for such a large range of particle sizes, from /xm to km. The values 
for smaller bodies on the order of a few kilograms are mostly determined from laboratory experiments, while the values 
for larger bodies are determined from under-surface ex plosions, observations of large asteroids and with experiments 
done under very high pressure (jHolsapple et al.l [20021 ). However, material strength varies greatly as a function of 
material type, object size, surface type and the number of shattering events an object has gone through over its 
lifetime. An object that has gone through many collisions in its lifetime, but still remains in one piece (descriptively 
called a "rubble pile") can endure harder collisions, which can actually be absorbed and help to compact the object, 
rather than dispersing it into smaller particles. This may seem like an important parameter only for larger objects; 
however, the evolution of larger objects significantly influences the evolution of smaller particles, and thus is important 
in our study. We also lack experiments done with targets and impactors cooled down to space temperatures of 100-150 
K, where one would assume that objects get more brittle and easier to shatter. 

Experiments clearly show that Q* is a function of the target mass M, meaning that d i fferen t mass targets will get 
shattered (with a 0.5M largest fragment) by different specific energies. iHolsapple et"all (|200l reviews experimental 
and theoretical results on collisions and strength curves. A common result for all of them is a minimum in the strength 
curve for bodies around 0.3 km in radius, where planetesimals are easiest to disperse (the number of cavities and 
cracks weakening the bodies increases, while self-gravitation is not dominant yet). As a result, there is a bump in 
the size distribution of minor planets in the solar system around this size. Smooth particle hydrodynamic (SPH) 
models give the Q* strength curve for larger bodies, while experiments help to anchor the curve for smaller rocks on 
the scale of a few cm in radius. It is still not clear whether the power-law shape of the curve can be extrapolated 
down to micron size particles, where experiments cannot be carried out. To study the collisional evolution of the 
smallest particles, the exact value of the strength curve must be known. In the absence of any models/experiments 
currently at those sizes, the best that can be done is a simple extrapolation of the strength curve to those regimes. 
IStewart fc Leinhardti (|2009[ ) introduce a velocity-dependent tensile strength curve, that is defined by variables such 
that it removes ambiguities over material density and projectile-to-target mass ratio. Their tensile strength curve is 
ideal for low-velocity (1-300 m s~^) collisions, such as those found during planet formation or at large radii debris 
disks. However, their universal relationship does not hold for conditions that strongly depart from the catastrophic 
disruption regime (Qimp < O.IQ*). 

In our models we use the lBenz fc Asphaugl ()1999l ) dispersion strength curve. It is derived from SPH models, represents 
a reasonable average of all previous strength curves, and is impact angle averaged. This curve can be written as (all 
units are in SI) 

J gerg-i kg"^ , (Al) 



Q*{a) = 10- 



S 



a 



1 cm 



Gp 



a 



1 cm 
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where the fiducial values in the equation are given in Table [T] 

B. MASS CONSERVATION OF THE MODEL 

A crucial test of any collisional code is for it to conserve the initial total mass of the system. Since particles are 
removed at the low mass end, this behavior can be complicated to verify. However, a system can only maintain its 
total mass numerically, if its collisional equations are mass conserving analytically. Here, we prove that our collisional 
equation is mass conserving. 

The collisional equation can be written as 

dn{m) 
dt 



dm' n(m)n{'m')a{m, m') 

J poo 

dfi / dMn{fi)n{M)a{fi, M)S [X{fi, M) ~ m] 
J II 

dfi dMn{fi)n{M)a{fi,M)R{m;fi,M) , (Bl) 
J II 

where R{m; fj,, M) is the redistribution function to mass m from /i + M collisions, such that 

/•oo 

/ dmR{m;n,M)m^ n + M ~ X{fi,M) , (B2) 
Jo 

and S is the Kronecker function. Multiplying Equation (jBl[) by m and integrating over dm gives 

dM dn{m)m T ^ ' r W ' W M 

— — = — — I dm I dm n[m)n(m )a(m,m jm 

dt dt Jq Jq 

l*OC /"OO /"CO 

d^ / dMn{n)n{M)a{n, M) / dm5[X{fj,, M) - m]m 
J n Jo 

dn dMn{fi)n{M)a{fi,M) dmR{m; fi, M)m , (B3) 
Jo J II Jo 

where 



dmS[X{^i,M) -m]m = X{^i,M) , (B4) 



resulting in 



dM 
'dT 



dm J dm'n{m)n{m')a{m,m')m + j d^ / dMn{fi)n{M)a{^, M){fi + M) (B5) 
The first integral can be separated into two sections as 

) />oo 

dm / dm'n{m)n{m')a{m,m')m= / dm / dm'n{m)n{m')a{m,m')m 



dm / dm'n{m)n{m')a{m,m')m . (B6) 



Since a{fi^M) is a symmetric function, we can swap the limits of integration for m and m' in the second integral of 
Equation (IB6p and, after making a change of variables of m = /i and m' = M in the first and m' = /x and m = M in 
the second integral, the full equation becomes 



d/i / dMn{^l)n{M)[a{^l,M)^l + a{M,n)M -a{n,M){p + M)] . (B7) 



dM 
'dT 

Since the collisional cross section is completely symmetric, the integral itself becomes zero, thus proving that our 
equation is mass conserving. 

C. NUMERICAL EVALUATION OF THE MODEL AND VERIFICATION TESTS 

The integro-differential equation presented in §3 must be integrated over 40 orders of magnitude in mass space, 
contains a double integral whose errors can easily increase if not evaluated carefully, and is bundled in a differential 
equation that evolves the number densities of dust grains and boulders within the same step. These characteristics 
demand attention in its numerical evaluation. In the following subsections we explain the numerical methods used 
to evaluate each integral and the ordinary differential equation (ODE). We also present verification and convergence 
tests for our code, which explain why such precisions are really necessary. 
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a{n) a(n) a(^) 

1 |j.m 1 mm 1 m 1 km 1 |j.m 1 mm 1 m 1 km 1 |im 1 mm 1m 1 km 




logio(^i)(kg) Iog1o{^i)(kg) Iog1o{^l)(kg) 

Figure 14. The valu es of G(m, m') as a function of the colliding masses. The thick contour is for G{m, m') = 100, which is roughly equal 
to the r value used in IDohnanvil II1969I ). The panels give the contours as a function of collision velocities. The coUisional velocities of 0.5, 
1.0, and 1.5 km correspond to debris ring radii of 100, 25, and 10 AU around an A spectral type star, respectively. The G{m,m') 
parameter is strongly dependent on the coUisional velocity. 



C.l. Taylor series expansion of Ti 

First, we expand equation (|13p to use a Taylor series when m' <C m and m' < fixim). For this, we rewrite M in 
terms of m and m' as 

M ==m + m'G{m, m') , (CI) 

where m'G{m,m') equals the cratered mass, m is the largest X{AI,m') particle created and G(m, m') can be f ound 
by root finding algorithms. As written, G{m,m') can be related to the F parameter used by ^Dohn anvil (|1969( ). for 
which he used a constant value of 130 for 5 km s"""^ collisions. We plot the value of G(m, m!) as a function of fj, and 
M in Figure with the thicker solid line giving the c ontour of G{m, m ') = 100. This contour lies at sizes reasonable 
for experiments in laboratory conditions, which is whv IDohnanvil ()1969D used a value close to it. The positions of the 
contours are a strong function of the interaction velocities. The m' < /Ltx(w) integrand can be written as 

I{m,m') — f{m')w{m')(Tjij{m')m' ^a{m)^i f{m)w{m)aw{m)(l + Z] 



-f{M)w{M)(j^{M) [l + G(m, m')Z^] " |z + [l + G(m, m')Z^] } j (C2) 



where Z ~ a{m') / a{m) and /(m) and f{m') are dimensionless number densities that can be expressed as 

We rewrite this integrand as 

/(m, m') — J{m')w{m')aw {m')m' ^ a{m)'^ f {m)w{m)aw (m) ^(1 + Z)'^ 

fiMMMMM) ^ ^^^^ -'{Z+[1 + Gim, m')Z^] ^ V] . (C4) 

The Taylor series for the components are 

{1 + Zf = 1 + 2Z + Z^ 

= Ti (C5) 

(C6) 



and 
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[l + G{m,m')Z^]-ri Z+ ^1 + G{ni,m')Z3 



-1 + 2Z + + 

l2G{m,m') 
^ [ 3 



2r]G{m, m ) 



r]G{m, m') 



Z' 



Z^ -r]G{m,m')Z^ 



(C7) 



Both f(M)/f{m) and w{M)/'w{m) are close to 1, while cr^ (M) /(Tm (to) deviates from 1 as m approaches TOmax- In 
those cases, the ratio can be expressed as 



Exp(^ 



e 



(M-to) = 1-P 

Af=m Q [1 - Exp ^ 



Uw {m) dM 
since we know that M — m — m'G{m, to'). We write this ratio as 

= 1 - . 

The integrand then takes the form 

I{m,m') — f{m')w{m')aw{fn')ni' ^ a{m)'^ f{m)w{m)a^^(m) x 
Rearranging it gives us 

1(171,111,') — /(to')w(to')(Tu,(to')to' ''a(TO)'^/(TO)u)(TO)cru,(TO) X •|ti 

When 



)] 



m'G{m,m') , (C8) 



w[m)j[m) 



(C9) 



(CIO) 



l-(l-i^) 



w{M)f{M) 
w{m)f{m) 



-X.^^^^(l-F) 



w{m) f {m) 



^ _ w{M)f{M) ^ ^^_9 ^ 
w{m)f{ni) ' 

we use the approximate formula 

I{m,m') — f{m')w{m')ayj{m')m' ^ a{ra)^ f{M)w{M)ayj{ra) x 



XiF~X2(l-F) 



(Cll) 
(C12) 

(C13) 



We use the Taylor series of the components to write the integrand below the limit of Z < 10 ^ (i.e., m! jm < 10 ^) 
This means that our full integral for the first term (Tj) takes the final form 



<ifiim,t) 
dt 



fj.x (m) 



-VtxG 



dm' I + J dm' f{m' ,t){m') J{m,t) {a{m) + a{m')Y 

fj.x (m) 



(C14) 



where 



'' f{m')w{m')aw{m')m' a{m)'^ f {m)w{m)(T^{m) x 



Ti 



1 _ (1 _ ^) !£LMiZiMi 

V / 10(771) J (m) 



_^ w(M)f{M) 
^ w(m)f(m) 



(l-F)} if to' < TO X 10-y & 1 



-9 T _ w(M)f(M) 
w{ra) f (rn) 



1= < 



f {m')w{m')aw{m')m' ^a{m)^ f{M)w{M)aw{m)x 

f{m', t)(TO')-'' [/(to, t) {a{m) + a{m')f - 

f{M,t){M.)-\a{M)+a{m')f 



if to' < TO x 10-9 & 1 - '"^.^'^/.^^'^ < 10-9 

w{rn)f{m) 



if TO X 10 ^ < m' < ^xifn) 



In Figure[Tni we show the mass of the X fragments created when particles of mass /i and M collide. The to = X{ii, M) 
regions are well defined in our single coUisional velocity case. When using collisional velocity that depends on particle 
size, more than one nxim) boundary may exist. 
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Figure 15. The largest X fragment produced by collisions between particles fi and M as a function of collision velocities. The collisional 
velocities of 0.5, 1.0, and 1.5 km correspond to debris ring radii of 100, 25, and 10 AU around an A spectral type star, respectively. 



C.2. Verification of the numerical precision of Tj 

To verify the precision of our integrator we set up an equation that is similar to Tj in behavior and that has an analytic 
solution and compare values given by our code to it. The integral we evaluate both analytically and numerically with 
our code is 



TVi{m) = 



{/ 1 i\2 r 1 ii2/TO + m'r\ 



+ 



(C15) 



where we have removed all the constants and the dimensionless number densities. We have also replaced G{m',m) 
with a constant F (which can be related to the F parameter used by Dohnanyi 1969) to enable an analytic solution. 
To verify our algorithm for the evaluation of this term we set the initial particle distribution to a power-law {rj—11/6). 
The integration boundary can be evaluated as 



if F < 1 
if F > 1 



The first integral is (setting Z = ^) 
2 




15Z3 + 10Z3 + 3 _ (1 + rz) 



10Z3(1 + 2ZF) 3(l + 6Zr) 3Z3(5 + 6Zr) 



{l + ZTY 



[i + zvy 



1 + ZT 



(C16) 



(C17) 



When ^ = Z < 10~^, the equation does not describe the analytic result correctly as the first three components 
completely cancel the last three components numerically; however, analytically the result is not zero and this non-zero 
component gets multiplied by a large number. This is the same catastrophic cancellation that affects the numerical 
evaluation of Ti(m). To overcome this and correctly represent the analytic result of the integral in such cases, we 
rewrote this to a Taylor series as well. The first three components cancel, and we are left with 



^1/6 



35 



TZ + IbTZ^'^ 



11 

Y 



65 
24 



jj^^5/3 ^2 j^2^7/3 



25, 

T 



85 
24' 



665 
432 



-r^z^ 



The second integrand has a much simpler anti-derivative 

6 m^/-^ 



5 



iV3 

7T72 



(C18) 



(C19) 



In Figure I16[ we show the computational error as a function of mass m, the F constant and the number of grid 
points used (neighboring grid point mass ratio). In the actual model F is not a constant, but equal to the variable 
G{m, to') which, as shown in Figure[T31 varies from 10"-* to lO^. Fi gure shows that the errors do not improve much 
past N=1000 {S = 1.1) and that, in general, errors are smaller for large values of F. This shows that errors originating 
from Ti are most likely to affect the smallest particles in the system. We expect the errors to actually be completely 
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Figure 16. The error in the integration of Ti as a function of the mass m and the value of the F constant for neighboring mass grid ratios 
of <5 =1.64, 1.18, 1.10 and 1.05. 



a{n) a(^l) a(^) 
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Figure 17. Iso-size contours for the produced Y fragments as a function of the colliding body sizes and interaction velocities. Fragments of 
size a(m) will be produced within regions where a{m) < a{Y). The largest fragments produced are not heavily dependent on the interaction 
velocities. The coUisional velocities of 0.5, 1.0 and 1.5 km s^^ correspond to debris ring radii of 100, 25 and 10 AU around an A spectral 
type star, respectively. 



symmetric, with the highest masses showing the same quantitative errors as the lowest masses near the boundary. 
Offsets are due to the fact that our analytic model includes targets larger than mmax- The maximum error of 10~^ is 
not ideal, but acceptable. When running our code, we use N=1000 grid points, which corresponds to & 6 — 1.1. 

C.3. Numerical evaluation of Tjj 

As a double integral, the second term, Tu, poses a larger challenge to achieve acceptable precision. A collision 
between masses fi and M will be able to produce a mass m in the redistribution power-law, if to < Y{^, M). In Figure 
[T7|we plot the iso-K contours in the /i vs. M phase space, which shows that integrating between exact boundaries is 
difficult for T//, especially if the coUisional velocity is not a constant but a function of the particle mass. 

As a first step, we determine which m masses can be produced by the grid points (/x, M) and their neighbors. For 
a grid point to be able to produce a particle of mass to, M) has to be larger than to. We determine the limiting 
mass that is produced by each grid point and all of their neighbors as well. Up to that min(/i, M) value, all to masses 
are produced with the full weight of the grid point. Between min(/Lt, M) and max(/x, M) - which is the largest to still 
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produced by the (/Lt, M) grid point itself - we analyze the areas divided into quadrants, and assign integration weights 
appropriately. A simple plot is shown in Figure [18] to explain the weights given to each grid point. 

These minimum and maximum masses are usually at most 2-3 grid points apart. This means that on average 2-3 
numbers have to be stored for all {fi; M) grid points, as below min(/i, M) all m masses have the same weight. The 
final integration speed can be increased by factors of 5 as the integration loops can be run in non-redundant ways. 

C.4. Convergence tests 




logio(M) 

Figure 18. Description of the integration method used for Tjj. The blue line represents the boundary, within which collisions are able 
to produce a certain mass m in the redistribution power-law. Resolution elements capable of producing m on their full area are colored 
red. Boundary resolution elements (i.e., fj, = M or M = rrimax) will be able to produce m, on a "half area, colored blue. The tip of 
the distribution (green) will be able to produce on an eighth of its full area. Partial quarter contributions are given by the yellow areas. 
Increasing the number of grid points used obviously increases not just the precision but the area used for the integration as well. 
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Figure 19. Convergence test results of our code. The left panel gives the results for Ti, while the right panel for T//. Wc also plot a A^^^ 
curve with a dashed line, which is the effective accuracy of the trapezoid integration method. The accuracy of the first term follows this 
trend, however, that of the second term is shallower, due to the resolution dependent integration limits. 

We run convergence tests on our code for both terms to find the least number of grid points we can use and still 
keep an acceptable numerical accuracy. We calculate convergence to the value given using 4000 grid points. Our 
convergence plot for Ti in Figure [19] shows that, for such a large dynamic range in masses, one needs a neighboring 
grid mass ratio of at most <5 = 1.1 to reach a relative error of 10~^. 

Our convergence test for Tjj does not reach the same level of accuracy, as at (5 « 1.1 we reach a relative error of 
10~^ only. However, this error is driven by the resolution dependent integration limits and not the method itself. As 
such the number of particles added by Tjj will always be underestimated by a small amount. 
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C.5. The ODE solver 

Previous work (i.e., iThebault et"al] 120031: iThebault fc Augereaul [20071 : IKrivov et all 120001: ILohne et al.1l2008[ ) used 
only a first order Eulerian algorithm to solve the difi^erential equation. We are using a 4"^ order Runge-Kutta algorithm 
(RK4) . To verify the ODE solver we simply evolve the Poynting- Robertson drag term, whose analytic solution is 

7i(m,i) = n(TO,t = 0)exp ( —] . (C20) 

In the following, we verified the accuracy of the ODE integrator by setting f3 = 0.100 for the particles and using the 
solar system timescale of 400 years. 

Using the results from the code, we define the ratio of the particle density at some time t to the particle density at 
time zero, i.e., i?code = fi'm'Tt)/f(rn,t = 0), and compare it to the analytic result, i?exp = exp(— i/rpRD)- We then 
compute the fractional error between the numerical and the analytic solution. The left panel of Figure [20] shows the 
fractional error as a function of the time step At, evaluated at a time roughly t « rp^D (4000 years) for both the RK4 
and Euler method. 

For this particular set-up, the optimal time step is ^ 4 yr for the RK4. However, the optimal time step depends on 
the time t at which the fractional difference is evaluated because of the accumulation of round-off errors. Finally, the 
right panel of Figure [20] shows the amount of CPU time taken to calculate an RK4 step as a function of the number 
of mass grid points used on a Mac Pro4 with 2 2.26 GHz Quad Core Intel Xeon processors. 
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Figure 20. Left Panel: The fractional difference between the numerical and analytical results at a time roughly equal to 4000 yr as a 
function of the time step used. We show errors for both an Euler ODE solver and for our RK4 algorithm. We also plot a and a t"* 
curve to guide the eye. Right Panel: The amount of processor time needed by our code to complete an RK4 step as a function of the 
number of mass grid points used. 
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